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A catchment scale numerical model is developed based on the three-dimensional transient Richards 
equation describing fluid flow in variably saturated porous media. The model is designed to take 
advantage of digital elevation data bases and of information extracted from these data bases by 
topographic analysis. The practical application of the model is demonstrated in simulations of a small 
subcatchment of the Konza Prairie reserve near Manhattan, Kansas. In a preliminary investigation of 
computational issues related to model resolution, we obtain satisfactory numerical results using large 
aspect ratios, suggesting that horizontal grid dimensions may not be unreasonably constrained by the 
typically much smaller vertical length scale of a catchment and by vertical discretization requirements. 
Additional tests are needed to examine the effects of numerical constraints and parameter heteroge- 
neity in determining acceptable grid aspect ratios. In other simulations we attempt to match the 
observed streamflow response of the catchment, and we point out the small contribution of the 
streamflow component to the overall water balance of the catchment. 


1. Introduction 

The grid resolution required to obtain acceptable numeri- 
cal solutions will be an important controlling factor in the 
computational feasibility of running large-scale catchment 
simulations. Discretization constraints can arise from phys- 
ical or numerical considerations. Time steps can be numer- 
ically constrained in order to satisfy convergence, stability, 
or accuracy requirements. A physical constraint would 
typically be one which is imposed in order to capture the 
dynamics of a process of interest, for instance, a time step of 
the order of seconds or minutes if one is interested in the 
timing and magnitude of surface saturation and runoff re- 
sponses during a heavy rainstorm. In the spatial domain one 
of the important constraints for a catchment scale subsurface 
model is connected to the aspect ratio of the numerical grid, 
which we define as the ratio of the horizontal mesh size A x , 
Ay to the vertical discretization A z. The horizontal extent of 
a large catchment will typically be much greater than its 
vertical length scale (large surface area and comparatively 
thin soil zone), and often a very fine vertical resolution, 
especially near the surface, is needed to accurately simulate 
infiltration and evaporation processes. Numerical catchment 
simulations are therefore computationally feasible so long as 
we can use grids with a large aspect ratio. If we are 
constrained to adopt smaller aspect ratios (decreasing Ar, 
Ay) in order to overcome numerical difficulties, then the size 
of the problem (number of degrees of freedom) can quickly 
exceed the capacity of available computers. 

In this paper we describe a physically based three- 
dimensional finite element model for the simulation of sub- 
surface hydrologic processes at the subcatchment and catch- 
ment scales, and we apply the model to a subcatchment of 
the Konza Prairie Research Natural Area near Manhattan, 
Kansas. The practical application of our model to actual 
catchments is demonstrated, and we investigate computa- 
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tional issues concerning the effects of model resolution 
(discretization, aggregation, and convergence constraints) 
on large-scale simulation of hydrologic processes. 

One of the overriding problems in hydrology is the under- 
standing of responses over the range of scales O(10 “ 1 ) to 
O(10 4 ) km 2 [Wood et al., 1988; Goodrich and Woolhiser , 
1991]. The low end of this range is roughly the size of a small 
subcatchment and marks a transition from point and hill- 
slope scales, at which physically based hydrologic models 
are more easily tested and better understood, to catchment 
or basin scales. The high end corresponds to the horizontal 
grid scale used in general circulation models for global 
climate simulations, and a better understanding of hydro- 
logic processes at this scale is required to improve the land 
surface boundary conditions for these models. The parame- 
terization of hydrologic processes at large scales is made 
difficult by the high degree of nonlinearity and variability in 
catchment parameters and inputs, and thus conceptual or 
idealized models are often used at these scales. Physically 
based analytical or numerical models can be used to study 
the validity of simplifying assumptions in conceptual mod- 
els. Some examples can be found in the works by Reeves 
and Miller [1975] (time compression approximation for par- 
titioning rainfall into runoff and infiltration), Broadbridge 
and White [1987] (time to ponding), Gan and Burges [1990a, 
b] (catchment rainfall-runoff models), Shamsai and 
Narasimhan [1991] (Dupuit-Forchheimer assumption under 
seepage face conditions), Sloan and Moore [1984] and Stag- 
nitti et al. [1986] (subsurface flow models), Wilcox et al. 
[1990] (runoff prediction models), and Troch et al. [1993] 
(catchment scale water balance models). 

Other studies using physically based hiilslope and catch- 
ment scale models include the early work of Freeze [1971, 
1972a, b ], who used a three-dimensional finite difference 
variably saturated flow model coupled with a one- 
dimensional channel flow model to reveal the importance of 
subsurface flow processes and parameter variability on 
watershed runoff response. Smith and Hebbert [1983] sim- 
plified the model used by Freeze and applied it to an 
experimental hiilslope in Western Australia to investigate 
the effects of rainfall, soil properties, and hiilslope geometry 
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Fig. I. Map showing location of the Konza Prairie Research 
Natural Area, the Kings Creek catchment, and the ID subcatch- 
ment. 


on runoff. Loague and Freeze [1985] and Loague [1990] 
compared the performance of a simple linear regression 
model, a unit hydrograph model, and a quasi-physically 
based model in simulating rainfall-runoff response on small 
catchments. Beven [1977], Bathurst [1986], and Govindaraju 
and Kavvas [1991] coupled one- or two-dimensional subsur- 
face flow models to physically based models of overland flow 
and channel flow. Binley et al. [1989a, b] used a three- 
dimensional model of variably saturated flow to explore the 
effects of spatially variable hydraulic conductivity and the 
validity of using an equivalent or effective conductivity 
value. 

The formulation of the numerical model presented here is 
consistent with one of the long-term objectives of our work, 
which is to simulate a large catchment such as the Kings 
Creek catchment shown in Figure 1, containing many dis- 


tinct subcatchment units and a complex stream network. In 
this paper we describe simulations of the ID subcatchment 
(Figures 1 and 2) located near the Kings Creek catchment 
(note that “ID” is a site name and bears no relation to the 
dimensionality of the subcatchment). The primary input to 
the model is digital topographic data obtained from the U.S. 
Geological Survey. The original data are in regular grid form, 
with a resolution of 30 x 30 m, but can be readily interpo- 
lated or extrapolated to finer or coarser discretizations. An 
automated extraction algorithm [Band, 1986] is applied to 
the topographic data to produce a mapping of stream chan- 
nels, catchment boundaries, and subcatchment partitions. 
Field observations and remotely sensed data are used for 
parameterization and calibration of the model. Atmospheric 
inputs to the model are specified as rainfall and evaporation 
boundary conditions, and the switching to and from atmo- 
sphere- and soil-controlled surface inputs is handled auto- 
matically. 

The catchment simulation model is based on the three- 
dimensional Richards equation describing fluid flow in vari- 
ably saturated porous media. Of the many numerical issues 
associated with large-scale three-dimensional simulations, 
some which are specific to catchment subsurface flow mod- 
eling include, first, the nonlinearity of Richards’s equation, 
and, at the interface between the saturated and unsaturated 
zones, possible discontinuities in the nonlinear coefficients 
and a change in type of the governing partial differential 
equation (parabolic to elliptic). The nonlinear system inte- 
grals require numerical evaluation or some other approxima- 
tion technique, introducing additional error in the model. 
Moreover convergence of iterative schemes used to solve 
the nonlinear equation cannot be guaranteed, and the rate of 
convergence will depend upon many factors (such as the 
time step, which directly affects the quality of the initial 
solution estimate used in an iterative procedure). A second 
issue is naturally occurring spatial and temporal variability in 
soils, topography, vegetation, rainfall, and evaporation, re- 
quiring complex boundary conditions and a high degree of 
heterogeneity in the model parameterization. This can pro- 
duce ill-conditioned system matrices and adversely affect the 
convergence behavior of linear and nonlinear iterative solv- 
ers. A third issue is the irregular geometry of catchments, 
resulting in sparse system matrices which are not regularly 
structured. A final issue to note is the large horizontal extent 
of a catchment compared to its vertical range and vertical 
discretization requirements, producing distorted elements. 

A series of test simulations on three small hypothetical 
catchments was conducted to evaluate the performance of 
the numerical code [ Paniconi , 1991]. The simulations in- 
volved alternating episodes of rainfall and evaporation, and 
generated significant amounts of discharge, infiltration, and 
saturation excess runoff. The tests included a comparison of 
lumped and distributed mass matrix versions of the code and 
a comparison of direct and iterative solvers for the linearized 
system of equations in the model. Various storage schemes 
for the system matrices and for the Jacobian coordinate 
transformation components were also examined. 

In the next section we introduce the catchment simulation 
model and describe in detail the generation of the numerical 
grid and the representation of various hydrologic processes. 
This is followed by a description of the Konza Prairie ID 
catchment and observation data. The simulation model is 
then applied to a 17-day rainfall-interstorm sequence where 
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Fig. 2. Elevation image of subcatchment I D- 10 ( 10 x 10 m grid resolution) showing stream (solid pixels) and location 
of three surface nodes selected for vertical profile output (shaded pixels). 


we discuss some calibration results. We use detailed and 
averaged rainfall rates for the 17-day simulation to illustrate 
temporal aggregation effects. Model resolution effects are 
discussed in more detail in section 4, where we use a 9-day 
interstorm simulation to examine aspect ratio and conver- 
gence constraints. 

2. Description of the Model 
2.1. Assumptions and Limitations 

Basing our numerical model for flow in variably saturated 
porous media on Darcy’s law and Richards’s equation, we 
adopt the usual set of assumptions: Flow is laminar and 
isothermal, inertial forces and chemical gradients are ne- 
glected, and the air phase is continuous and at atmospheric 
pressure [e.g., Freeze and Cherry , 1979; Hillel, 1980a; 
Sposito , 1986]. In addition, we do not account for hysteresis, 
we assume that the porous medium is isotropic, and we 
consider only flow within the “soil matrix,’’ neglecting flow 
through “macropores.” Anisotropy can be easily incorpo- 
rated with a generalization of the hydraulic conductivity 
term in the model equations. Since our model treats both 
precipitation and evaporation inputs, it will be important to 
consider hysteresis effects in future versions of the code. 
Mualem [1974], Parlange [1976], and Kooi and Parker [1987] 
discuss some conceptual models of hysteresis for soil mois- 
ture characteristic equations. 

Whereas anisotropy and hysteresis can be readily incor- 
porated into our catchment simulation model, treating non- 


isothermal effects and macropore flow would require signif- 
icant extensions or generalizations of the model. Several 
recent studies suggest that macropore flow (also described as 
bypass flow, channeling, pipe flow, or preferential flow) may 
contribute significantly to the transport of water through 
hillslopes and catchments [e.g., McDonnell , 1990; Pearce , 
1990]. The term macropore generally refers to continuous 
pore structures which can exhibit nonequilibrium channeling 
flow, and it may not be appropriate to model this type of flow 
based on Darcy’s law [Beven and Germann , 1982]. Two- 
domain models have been proposed for describing the flow 
and interactions in a soil matrix/macropore system, and 
kinematic wave models have also been introduced [e.g., 
Germann , 1990]. 

Subsurface heat transport and temperature-induced soil 
moisture flow are important not only for long-term simula- 
tions which need to account for seasonal changes in temper- 
ature, but also for short-term simulations when one consid- 
ers the effects of diurnal fluctuations in solar and 
atmospheric radiation on near-surface soil moisture pro- 
cesses such as evaporation. Milly [1982] presents a physi- 
cally based one-dimensional coupled moisture and heat flow 
model based on studies of nonisothermal flow in porous 
media by Philip and de Vries [1957] and de Vries [1958]. The 
model treats moisture flow in both the liquid and vapor 
phases and latent and sensible heat transport by conduction 
and advection. 

The limitations described above pertain to the physics of 
flow processs within a porous medium. In a realistic basin 
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scale model one must consider not only the coupling of 
processes in the saturated and unsaturated zones (soil matrix 
and macropores, air and liquid phases, heat and moisture 
flows), but also the coupling between subsurface, surface, 
and atmospheric hydrologic processes. Surface processes 
which affect and are affected by subsurface moisture and 
energy states include overland flow and streamflow [ Eagle- 
son , 1970; Freeze , 1974] and vegetation growth and transpi- 
ration [e.g., Federer , 1979]. The coupling with atmospheric 
processes takes us into the realm of large-scale water and 
energy balance models [e.g., Eagleson, 1978, 1979] which 
can become a component of general circulation models. 
Avissar and Verstraete [1990] discuss some of the difficulties 
involved in representing small-scale surface processes (such 
as albedo, stomatal response, momentum transfer, surface 
roughness, and soil moisture flow) in large-scale atmospheric 
models. 

2.2. Governing Equations and Numerical Procedures 

The three-dimensional Richards equation with pressure 
head ^ as the dependent variable can be written as 

dtff 

SW-=V-[K s K r WV(<l' + z)] ( 1 ) 

ot 

where / is time, z is the vertical coordinate, positive upward, 
and the hydraulic conductivity K is expressed as a product 
of the conductivity at saturation, K s> and the relative 
conductivity, K r . We use an extension of the van Genuchten 
characteristic equations [van Genuchten and Nielsen , 1985] 
to describe the nonlinear dependencies of volumetric mois- 
ture content 0, specific moisture capacity 5, and relative 
hydraulic conductivity K r on pressure head [Paniconi et aL, 
1991]: 


To solve (1) numerically we use a finite element Galerkin 
discretization in space and a finite difference discretization 
of the time derivative term [e.g., Ames, 1977; Huyakorn and 
Pinder , 1983]. The problem domain is discretized into M 
hexahedral elements and an approximating function is intro- 
duced: 


M 

ik(x, y, z, /) - Hx, y> z, t) = ^ y> z * 

e - I 


In local coordinate space (f, tj, £) the approximating function 
for each element ( e ) is written as 

8 

V, (, I) = 2 V, ()<pl e) «) 

i - 1 


= (N w (£, TJ, fl) r * w </) (7) 

where are undetermined nodal values of <// and N}*’ are 
trilinear Lagrange basis functions which define, in local 
coordinate space, a cubic element with eight nodes at (±1, 
± 1 , ±1). The basis functions have the general form N,(f, tj, 
0 = (1 ± f)(l ± tj)( 1 ± £)/8. The finite element 
formulation used is isoparametric, with the mapping from 
local (f, tj, 0 to global \x, y, z) coordinate space given by 

(*, y, z) = N >(Z' T »’ £)*/. 

8 8 

7V,(f, Tj, C)y h 2 *;(£. V, ()z, 

i - i /-i 



«(</-)= e r + (#, - e r )[i + p]~ m o 

( 

o(<f>) = e, + (o, - « r )[i + p 0 r m + - * 0 ) 

tl>> il> 0 


v dO (n- 



* d* |^,|"(1 + 1 

l/r < l/<0 

(3) 

dO 



W) = T7 =s, 

difi 

<t> ^ <£o 


K,(+) « (1 + |3)- 5m/2 t(l + /3) m - /3 m ] 2 

i/> < 0 

(4) 

K r (+) = 1 

<l> & 0 


where 0 r is the residual saturation, 6 S 

is the saturated 


moisture content, *f/ s is the capillary or air entry pressure 
head value, S s is the specific storage, m = 1 - Mn , 0 “ 
((/r/^) w , 00 - 0 (<Ao) = is a continuity 

parameter, and n can be interpreted as a pore size distribu- 
tion index. The exponential relationship 

K, = K s (z) = *s 9 exp [-AL - z)] ( 5 ) 

is used to model vertical heterogeneity of saturated hydrau- 
lic conductivity [Beven, 1982, 1984], where K s # is the satu- 
rated conductivity at the surface, /is a fitting parameter, and 
L is the elevation at the surface above the datum z = 0. 


where (x it y>, z*), / “ 1, * * * , 8 are the global coordinates 
of the eight comer nodes of element (e). 

The finite element spatial discretization yields the system 
of ordinary differential equations 


dV 

A (V)¥ + F(¥) — = q(0 - b(¥) 
dt 

(9) 

where ¥ is the vector of undetermined coefficients repre- 
senting the value of pressure head at each node. The system 
components for each element (e) are 

A (,) = JjJ K< e) K,(i (e) )[ Nj?(Nj‘J) T 


+ N^N^V+N^N^V] dCl (e> 

(10) 

b (,) = JJJ K< e) K r (ij> M ) d(l M 

(11) 

F (‘) = JJJ S(^ (e) )N M (N (e> ) T dO M 

(12) 

q(,) = f [ 9r«N (rt drtf 

J Jr„ 

(13) 
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where > is the specified Darcy flux on the natural bound* 
ary rj^. We use the notation N ?x to denote differentiation 
with respect to x, and N 7 to denote the transpose of N. A 
lumped form of the mass matrix F (e) is used: 

S(i/f M )N le) dCl ie) (14) 

where now N u) is taken to be a diagonal matrix rather than 
a vector. 

The system integrals are nonlinear and must be evaluated 
numerically. We use order 2 Gaussian quadrature with 
weights of 1.0 and Gauss points at (£, 17, £) = (±1/V3, 
± 1/V3, ± 1/V3). The saturated hydraulic conductivity K s 
is assumed constant over each element. 

To evaluate the system integrals in local coordinate space 
we need to compute the determinant of the Jacobian of the 
global to local coordinate transformation. For the integrals 
(10) and (11) we also need the inverse of the Jacobian to 
compute the derivatives of the basis functions (for instance, 
dNfdx = (dNild€)(d€!dx) + (dNi/drj)(d7)/dx) + (d N f 
The components of the Jacobian and inverse 
Jacobian will be spatially dependent due to the quadratic and 
cubic terms in the basis functions. The cost of computing 
and storing the Jacobian determinant and inverse Jacobian at 
each Gauss integration point for each element can be quite 
high for a three-dimensional simulation. For our particular 
model we could have taken advantage of special features of 
our catchment discretization to simplify the transformation 
of (8), thereby minimizing the storage and CPU expenses 
associated with the Jacobians. For instance, some of the 
spatial dependencies in (8) can be eliminated when using 
rectangular elements with the element edges aligned with the 
global coordinate axes. Furthermore, there are similarities in 
the Jacobian components from element to element which can 
be exploited if we use grid spacing in any of the coordinate 
directions. In the current implementation of the model we 
use the most general form of the mapping from local to global 
coordinates, as given by (8), without taking advantage of 
particular features of our grid geometry. 

The Crank-Nicolson finite difference scheme is used to 
discretize (9) in time. The resulting equation is 

* + t _ * 

A (X| t k + + U/2) + + (l/2)j 

A / 

= q(f* + ll/2) ) - b(^ i + <l/2 ’) (15) 

where superscript k represents time level and 'F* +<l/2> = 

Equation (15) is linearized using Picard iteration, which 
we can write as 

A * + (1/2), (m) + J_ Y k + * + l,(«f+ I) 

_ ij/* + U«)) - + h(m)j 

where superscript (m) denotes iteration level and 

f('I f * + S - + H/2)^^A + 1 + 

+ F* + < 1/2 > — - q* + »' 2 > + b* + itm = 0 (17) 

A / 


To solve the linear system of equations represented in (16) 
we use a conjugate gradient algorithm from ITPACK 
[Kincaid et al., 1982] with a symmetric successive overre- 
laxation preconditioner and a compact nonsymmetric stor- 
age scheme. The structure of the system matrix resulting 
from (16) is sparse and symmetric with a maximum of 27 (3 3 ) 
nonzero entries per row. For irregular catchment geometries 
such as shown in Figure 2 the bandwidth will not be constant 
over the matrix and can become quite large for some rows. 

2.3. Model Inputs and Representation 
of Hydrologic Processes 

The catchment simulation model comprises two programs: 
a grid generator which constructs the finite element mesh 
and initializes various parameters, and the actual simulation 
program which numerically solves the three-dimensional 
Richards equation over a specified time period for a given set 
of boundary and initial conditions. 

The design and structure of the mesh generation and 
simulation programs were motivated in large part by the 
availability of digital elevation models (DEMs) of topogra- 
phy from the U.S. Geological Survey. These data bases 
provide topographic information for extensive geographic 
regions at high resolution and in regular grid form (30 x 30 m 
pixels). An automated extraction algorithm [Band, 1986] is 
applied to the DEM data to produce a mapping of stream 
channels, ridges, and drainage basins and subbasins. Based 
on this mapping we can define the boundaries of a catch- 
ment, obtain elevation, stream, and distance-to-stream data, 
and subdivide the catchment into physically consistent sub- 
catchments. Figure 2 is a subcatchment image produced 
using some of the output from the extraction algorithm. The 
three shaded pixels in this image show the location of three 
surface nodes selected for detailed vertical profile output 
during the numerical simulation and the solid pixels outline 
the stream network. 

The finite element mesh generator was developed to take 
advantage of the regular grid structure of the digital elevation 
data and to use the information provided by the extraction 
algorithm. The mesh generation algorithm discretizes a 
catchment or subcatchment into hexahedral elements, num- 
bers and connects the nodes and elements, initializes the 
pointer arrays for storing the system matrices, and sets up 
the boundary and initial conditions for ensuring simulation. 
A compact storage scheme is used for the sparse system 
matrices. 

The simulation program takes the information from the 
grid generator, and for each time step of the simulation 
period it performs the iterations on the nonlinear equation, 
sets up and solves the linearized system of equations, 
calculates mass balance errors, and computes the hydro- 
graph contributions. The simulation program is an extension 
of the model by Binley et al. [1989a, b). 

The geometry of a catchment as defined by the DEM 
extraction algorithm is based on the location of naturally 
occurring ridges. The finite element method allows us to 
model irregular domains so there is no smoothing or trans- 
formation applied to redefine the catchment boundaries. The 
ridges are assumed to represent vertical walls which we 
consider to be impermeable lateral boundaries, and we also 
define the base of the catchment as a no-flux boundary. The 
only boundary conditions which need to be explicitly input 
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are those at the surface (precipitation and evaporation rates). 
This method of treating nonsurface boundaries is a simplifi- 
cation of the acutal physical processes, and its validity needs 
to be examined in more detail. In principle the model can 
handle any type of boundary condition, but there is often a 
lack of knowledge and data concerning flow processes 
across nonsurface boundaries. If we define the vertical 
extent of a catchment deep enough so as not to affect 
processes (infiltration, runoff, evaporation) occurring near 
the surface, or if we know the location of an underlying 
low-permeability layer, then an impermeable boundary con- 
dition at the base of the catchment is justified. Defining a 
deep catchment and analyzing, for example, the simulated 
flow behavior near the water table could be one method of 
determining more suitable boundary conditions when a 
no-flux condition is inappropriate at the catchment base. A 
similar method for assessing the validity of the zero-flux 
condition along lateral boundaries would be to simulate a 
catchment and examine the computed flow patterns across 
the ridge boundaries of each of the subcatchments which 
make up the catchment. 

The catchment topography is heterogeneous, with eleva- 
tion inputs obtained directly from the DEM data. The data 
are in regular grid form, and we retain this uniform grid 
structure in defining the x and y (horizontal) dimensions of 
each element of our numerical grid. The size of each DEM 
pixel is A x x Ay = 30 x 30 m, and it is an easy matter to 
interpolate these data to a finer grid or to aggregate it to a 
coarser resolution. In our simulations of the ID catchment 
we used linear interpolation to obtain a range of catchment 
discretizations from 30 x 30 m to 1 x 1 m. Although we use 
uniform grid spacing horizontally (Ax = Ay, constant), the 
vertical discretization A z can be variable. This allows us to 
define, for instance, thinner layers closest to the surface. 
The mesh generation program has the option of making the 
thickness of each layer uniform horizontally (this results in a 
staggered catchment base), or of making the base of the 
catchment flat (resulting in horizontally nonuniform layer 
thicknesses). 

To treat flow processes occurring on the catchment sur- 
face and in streams we need to couple Richards’s equation 
governing subsurface flow, the shallow water equations for 
channel flow, and a kinematic or shallow water model for 
overland flow [ Eagleson , 1970; Freeze, 1974]. In the current 
version of our model we use a simple linear transformation 
to distribute overland flow. Surface runoff generated at a 
point on the catchment is routed to the stream via a time 
delay determined from the overland flow velocity (currently 
assumed constant for the catchment) and the shortest dis- 
tance from the point to the stream. The location of the 
stream and values for the shortest overland paths to the 
stream (“delay distances”) are obtained from the DEM 
extraction algorithm. We note that this approach does not 
allow for downslope reinfiltration of overland runoff. This 
simple treatment of overland flow routing is a reasonable 
approximation when partial contributing areas are the main 
source of surface runoff. These saturated regions will typi- 
cally grow as an expanded stream network during a rainfall 
event [Dunne and Black , 1970] so that downslope reinfiltra- 
tion will not occur and the distances for overland travel will 
remain relatively short. Streams can be modeled as specified 
head boundaries (permanent streams), or we can assign 
stream nodes initially high levels of saturation (ephemeral 


streams). Channel flow is not considered in the model. 
Seepage faces and stream banks can be handled numerically 
[ Neuman , 1973; Cooley , 1983; Huyakorn et al. t 1986] but are 
somewhat complicated to discretize in the case where the 
stream is internal to the catchment rather than being situated 
along a lateral boundary. Wide internal streams require a 
dual stream bank configuration, with the location of stream 
banks and seepage faces automatically generated from anal- 
ysis of the DEM data. The current implementation of the 
model does not have this feature. 

In addition to heterogeneities in topography it is important 
to account for variability in the catchment soils and in the 
atmospheric inputs to the model. At present the model 
considers only spatial variability in saturated hydraulic con- 
ductivity K s and spatial and temporal variability in evapo- 
ration and precipitation rates. Other parameters, namely/, 
6 r , n, A,, and ^ min described below, are kept 

constant over the catchment. The extension to spatially 
variable representation of these parameters is straightfor- 
ward to implement in the model, although high levels of 
parameter heterogeneity may adversely affect numerical 
performance [Ababou et al. , 1989]. Values for the soil 
hydraulic parameters 6 r , 6 S , ip s , n , and S s can be obtained 
by fitting (2) and/or (4) to observed data or by using other 
information about the catchment soils. Saturated conductiv- 
ities are input for each node on the catchment surface, and 
relationship (5) is used to assign K s values vertically. The 
value of parameter / can be estimated by fitting observed K s 
data to (5), as shown later. By allowing spatial and temporal 
variability in precipitation and evaporation it is possible to 
simulate alternating periods of soil wetting and drying, and 
to have rainfall and evaporation occurring simultaneously 
over different portions of the catchment. For the simulations 
described in this paper we used spatially homogeneous 
rainfall and evaporation since the ID catchment is quite 
small. 

At any surface node the simulation program automatically 
switches from a specified flux (Neumann) to a constant head 
(Dirichlet) boundary condition when the node becomes 
saturated or its pressure head becomes smaller than the 
“air-dry” pressure head value [Hillel, 19805, p. 121]. 
The boundary condition switches back to a Neumann type 
when the magnitude of the flux across the soil surface 
(computed) exceeds the magnitude of the atmospheric (spec- 
ified) flux, or when the atmospheric event switches from 
rainfall to evaporation or evaporation to rainfall. For Di- 
richlet nodes the flux across the surface is computed by back 
solving (16) for q after having solved for the pressure heads. 
Surface boundary conditions are updated in this manner 
after each nonlinear iteration. 

The initial conditions required for a transient simulation 
are input as nodal pressure head values. The initial head 
distribution can be obtained by solving a steady state prob- 
lem, for example. Alternatively, we can generate the initial 
heads based on knowledge of the initial water table distribu- 
tion or initial soil saturation deficits. One method of calcu- 
lating water table depths or saturation deficits uses a topo- 
graphic index (computed from the digital elevation data) 
together with surface K s values, the exponential parameter 
/, and a base flow parameter [Sivapalan et a/., 1987]. A 
water table depth or soil moisture deficit can be converted 
into a vertical pressure head distribution using, for instance, 
a hydrostatic assumption. 
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The final group of input parameters to the simulation 
program are for dynamic time step control, back stepping, 
monitoring convergence of the nonlinear iterations, and 
controlling the generation of output for postprocessing. The 
specified convergence tolerance is to! , and the maximum 
number of nonlinear iterations allowed during a time step is 
maxit . In our simulations we used the convergence test 
|jt|r*+ 1.(« + D - ^ir*+t,(m)||^ < to f ^e simulation begins 
with a time step of Ar 0 and proceeds until time T m ax . The 
current time step is increased by a factor of Af mag (to a 
maximum of A/ max ) if convergence of the nonlinear system is 
achieved in fewer than maxit 2 iterations, whereas the time 
step is decreased by a factor of A/ rcd (to a minimum of A/ min ) 
if convergence required more than maxit 2 iterations. We 
back step with a reduced time step (factor A/ rcd , to a 
minimum of A/ min ) if convergence is not attained {maxit 
exceeded). 

2.4. Model Outputs 

Catchment simulation output is processed to produce 
plots, images, and summary statistics. Vertical profiles of 
pressure head and moisture content are plotted at various 
times to show, for instance, the water table response during 
the simulation. Shaded or color images of surface saturation 
can show the different mechanisms contributing to surface 
runoff on various portions of the catchment and the growth 
of partial contributing areas. Surface images of moisture 
content, pressure head, and flux values can also be easily 
produced, as can images along nonsurface cross sections of 
the catchment. 

Hydrograph plots of actual and potential catchment in- 
flows and outflows are produced. The potential inflow in a 
hydrograph plot consists of the precipitation (positive) and 
evaporation (negative) flux inputs supplied to the model. 
When the potential flux is positive, the difference between 
potential and actual soil inflow is the total runoff. Surface 
runoff is produced when the surface becomes saturated, 
either due to a rising water table (saturation excess mecha- 
nism) or to the infiltration capacity of the soil’s falling below 
the rainfall rate (infiltration excess mechanism) { Freeze , 
1974]. In both cases the boundary condition at the node on 
the surface where saturation occurs switches from a Neu- 
mann type (atmosphere-controlled inflow) to a Dirichlet type 
(soil-controlled inflow). Subsurface runoff in our model is 
only produced when the soil moisture flux becomes negative 
across the surface during a rainfall event, that is, when 
subsurface water exits the soil matrix from a saturated 
region on the surface. This type of subsurface runoff is 
sometimes called “return flow” [Dunne and Black , 1970]. 
The path of return flow to the stream is initially subsurface 
and becomes overland flow when it emerges at the surface. 
Without stream banks there is no way to treat seepage face 
subsurface flow. The stream discharges in a hydrograph plot 
are the surface and subsurface runoff components routed to 
the stream via the time delays computed from the overland 
flow velocity and shortest overland paths to the stream. 
There is no routing along the stream channel to the catch- 
ment outlet. 

The convergence behavior (for both the nonlinear iterative 
scheme and the linear solver) and mass balance errors 
(absolute and normalized) are plotted for each catchment 
simulation. These plots can help pinpoint trouble spots in a 


simulation, which often occur during transition periods 
between rainfall and evaporation, during heavy rainfall epi- 
sodes, or when the catchment or catchment surface become 
highly saturated. Absolute mass balance errors are com- 
puted at each time step as the difference between the change 
in moisture storage in the catchment soil and the net bound- 
ary influx (amount of water entering or leaving the catch- 
ment at the boundaries). The moisture storage is calculated 
by integrating the 0(i/4 (equation (2)) over each element and 
summing over all elements. In the current version of the 
model, nonzero boundary fluxes occur only at the catchment 
surface, so the net influx is obtained using the specified 
atmospheric rates (for Neumann nodes) or the back solved 
fluxes (for Dirichlet nodes). Normalized or relative mass 
balance errors are obtained by dividing the absolute mass 
balance errors by the net influx (and multiplying the result by 
100 for a percentage error). 

3. Simulation of the Konza Prairie 
ID Catchment 

We apply our three-dimensional numerical model to sim- 
ulate a subcatchment (catchment “ID”) of the Konza Prairie 
Research Natural Area in northeastern Kansas. Observation 
data collected during a 1987 field experiment are used to 
parameterize and calibrate the model. 

Whereas extensive streamflow, evaporation (actual rates), 
and rainfall data were available for the simulation periods of 
interest, very little soil data specific to the ID subcatchment 
were collected, and in some cases parameter estimates were 
made using generic soil characteristics for the region. Pa- 
rameterization of the model is therefore adequate for the 
atmospheric boundary conditions but less than satisfactory 
for the saturated conductivity distribution, soil zone depths, 
characteristic equations, and initial conditions. Due to lack 
of data some degree of calibration is needed even though the 
model is physically based. Moreover, the streamflow com- 
ponent of the ID catchment water balance is of much smaller 
magnitude than water losses due to evaporation, so that 
model calibration based on observed streamflow hydro- 
graphs alone is of limited use. For these reasons a compre- 
hensive and rigorous model assessment, according to proce- 
dures such as those described by James and Burges [1982], 
will not be our focus. Evaluation of possible model errors 
and a more detailed calibration should be conducted using 
potential evaporation measurements and soil moisture con- 
tent data in addition to streamflow observations. Data and 
measurement error would also need to be quantified and 
taken into account. 

3.1. Description and Discretization of the Catchment 

The Kings Creek catchment and ID subcatchment are in 
the Konza Prairie Research Natural Area near Manhattan, 
Kansas, in the northeastern part of the state (Figure 1). The 
Konza site, part of the Flint Hills Upland geologial region, is 
a 34.87 km 2 area of native tallgrass (blue stem) prairie and is 
one of 1 1 sites selected for the Long-Term Ecological 
Research (LTER) program of the National Science Founda- 
tion [ Bhowmik , 1987]. The Konza area has a temperate 
midcontinental climate with average annual precipitation of 
about 835 mm. Streamflow is ephemeral with high flows 
during the winter and late spring but dry in the summer and 
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TABLE 1. Discretization of ID Catchment 


Catchment 

Grid 

Resolution 

Number of 
Nodes 

Number of 
Elements 

Sparsity, 

% 

Storage, 
x 10 6 words 

ID-30 

30 x 30 m 

1,570 

1,064 

98.280 

0.148 

ID-15 

15 x 15 m 

5,795 

4,256 

99.534 

0.549 

ID-10 

10 x 10 m 

12,680 

9,576 

99.787 

1.20 

ID-06 

6 x 6 m 

34,430 

26,600 

99.922 

3.28 

ID-05 

5 x 5 m 

49,295 

38,304 

99.945 

4.70 

ID-03 

3 x 3 m 

135,355 

106,400 

99.980 

12.91 


late fall except during heavy rainstorms [Engman et at ., 
1989]. The ID subcatchment has a surface area of 0.24 km 2 
and is situated in the southeast comer just off the Kings 
Creek catchment. The U.S. Geological Survey digital eleva- 
tion data for the ID catchment contain 314 pixels at 30 x 30 
m resolution. 

Soil, streamflow, rainfall, and evaporation data at the 
Konza site were collected during summer 1987, from late 
May to mid-October, as part of an international field exper- 
iment (the First International Satellite Land Surface Clima- 
tology Project Field Experiment, or FIFE) to study land- 
atmospheric interactions for global climate modeling [Sellers 
et al., 1990]. Using the observation data we applied our 
catchment model to simulate the ID subcatchment for a 
9-day interstorm sequence (May 29 to June 6, 1987) and for 
a 17-day period of alternating rainfall and evaporation (June 
25 to July 11, 1987). 

To study aspect ratio constraints we interpolated the 30 x 
30 m elevation data for the ID subcatchment to obtain a 
range of discretizations from 30 x 30 m to 3 x 3 m, shown in 
Table 1. We label the 30 x 30 m discretized catchment 
ID-30, the 15 x 15 m ID-15, and so on. Figure 2 is an image 
of the ID catchment at a horizontal grid discretization of 
10 x 10 m. The vertical discretization of catchment ID was 
kept fixed at four layers for all simulations, and the number 
of nodes and elements in the three-dimensional finite ele- 
ment mesh using the different horizontal grid resolutions 
ranges from 0(1O 3 ) to 0(1O 5 ). We also interpolated the 
original elevation data to 2 x 2 m and 1 x 1 m resolution 
(yielding finite element grids of approximately 300,000 and 
1,200,000 nodes, respectively), but these grids could not be 
simulated as the amount of memory needed to run the code 
exceeded the 25 million word capacity (64-bit words) of the 
Pittsburgh Supercomputing Center Cray Y-MP computer. 
The degree of sparsity shown in Table 1 is calculated as 
100(1 - 27/jV), where N is the number of nodes and 27 is the 
maximum number of nonzero entries per row in the system 
matrices. 

3.2. Observation Data and Parameter Fitting 

The catchment and model parameter values used for the 
ID catchment simulations are given in Tables 2 and 3. The 
coordinates of three surface nodes designated for detailed 
vertical profile output are (x = 150 m, y = 360 m), (300, 360), 
and (450, 360) and correspond, from left to right, to the three 
shaded pixels shown in Figure 2. 

The values of «, /, and surface K s were obtained by 
least squares fits of (2) and (5) to observation data, as shown 
in Figure 3. The soils in part of the Kings Creek catchment 
are identified as Florence silty clay loams and Benfield silt 
loams, and we used moisture retention and saturated hy- 


draulic conductivity data from the Soil Conservation Service 
of the United States Department of Agriculture (T. Deme- 
triades-Shah et al. unpublished report, 1989) corresponding 
to Florence and Benfield soils at the FIFE sites. The values 
for residual moisture content 9 r and porosity 0* were taken 
from the estimates for a silty clay loam published by Rawls 
et at. [1982]. The curve and data in Figure 3 b suggest that 
even at extremely high suction values the soil moisture 
content does not fall below 20%. The measurements from 
the Soil Conservation Service shown in this plot were 
obtained from soil cores taken at depths of 5-10 cm and 
30-70 cm. From near-surface remotely sensed measure- 
ments of soil moisture made for the catchment, on the other 
hand, we do obtain moisture contents below 20% [Engman 
et al,, 1989]. Aside from measurement error, the discrepan- 
cies between the observations made by remote sensing and 
those obtained from soil cores may indicate differences in 
the characteristics of the soil near the surface and at depth, 
in which case these differences should be taken into account 
in the parameterization of the model. 

No measurements of relative hydraulic conductivity, spe- 
cific storage, or overland flow velocity were made for the ID 
catchment, and little information for estimating initial con- 
ditions or the horizontal distribution of saturated conduc- 
tivites was available. The K r (i/>) relationship was obtained 
from (4) using the parameter values from the 0(^) fit shown 
in Figure 3 b. We note that n = 1.176 is just outside the range 
1.25 < n < 6 suggested by van Genuchten and Nielsen [1985] 
for the validity of (4) and of the relationship m — 1 — Mn, 
so it will be important to validate the K r {\ji) relationship with 
experimental data. The horizontal distribution of saturated 
hydraulic conductivities was assumed to be spatially homo- 
geneous. For the 17-day precipitation-evaporation simula- 


TABLE 2. Soil and Grid Parameter Values for ID Catchment 
Simulations 


Parameter 

Value 

0 r 

0.04 

0 S 

0.471 

ift M > m 

-0.741 

n 

1.176 

S,, m _l 

0.001 

^min* ^ 

“-999.9” 

Surface elevation range, m 

[417.0, 442.0] 

A*, Ay, m 

3.0, 5, 6, 10, 15, 30.0 

(ID-03, ID-05, • *, ID-30) 

L, m 

1.0 

Az, m 


Top two layers 

0.2 

Bottom two layers 

0.3 

Overland flow velocity, m/h 

4.0 
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TABLE 3. Saturated Conductivity, Time Step, and Convergence Parameter Values for ID 

Catchment Simulations 


Parameter 

17-Day Precipitation- 
Evaporation Simulation 

9-Day 

Evaporation Simulation 

12-Hour 

Evaporation Simulation 

Surface K s , m/h 

0.2180 

0.0218 

0.0218 

/, m _l 

3.526 

3.526 

3.526 

Tmax. hours 

408.0 (17 days) 

216.0 (9 days) 

12.0 (0.5 days) 

Ar 0 , hours 

0.05 

0.10 

0.10 

A/ min, hours 

0.0005 

0.01 

0.0001 

A/ max * hours 

1.0*. 2.0t 

12.0 

3.0 

A/r*d 

0.5 

0.5 

0.5 

A/fymj 

1.2 

1.2 

1.2 

tol, m 

0.005 

0.05 

0.0005 

maxit 

8 

12 

12 

maxit 2 

5 

6 

6 


•Simulation using 15-min rainfall rates. 
tSimulation using daily-averaged rainfall rates. 


tion we assumed fairly wet antecedent moisture conditions we used a uniform value of 1.0 m in our simulations. We 

and tried various distributions of initial conditions. Varying discretized the catchment vertically into four layers with two 

the initial conditions, however, did not have as great an thin layers (0.2 m) near the surface and two thicker layers 

effect on the simulated catchment discharge rates as varia- (0.3 m) at the base of the catchment, 

tions in surface K s . For the 9-day evaporation simulation we Since catchment ID is quite small we assumed spatial 
assumed the catchment was completely saturated at the start homogeneity for the atmospheric inputs to the model. Pre- 

(based on high streamflows observed in the days just prior to cipitation data were collected at 32 rain gauges on the Kings 

the interstorm period, with a peak rate of approximately 26 Creek FIFE site, and we averaged the data from the two 

m 3 /h), and we used a vertically hydrostatic initial pressure gauges closest to the ID catchment. Evaporation measure- 
head distribution. We used values of 0.001 m" 1 for specific menls were made during four separate periods over the 

storage and 4.0 m/h for overland flow velocity, which will summer and represent actual rather than potential soil mois- 

require verification. The average soil depth of the ID catch- ture losses to the atmosphere. In the absence of potential 

ment has been estimated to be 0.76 m [ Blain , 1989], although evaporation rates it was necessary to set the air-dry pressure 



0.0 0.005 0.010 0.015 0.020 0.025 -150 -100 -50 0 


Saturated Hydraulic Conductivity (rnfaj Pressure Head (m) 

Fig. 3. Least squares fits (solid lines) of the (a) K s {z) and (b) 0(i/O relationships to observed data (crosses). 
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Time (day of year) 

Fig. 4. Observed precipitation data for the 17-day rainfall- 
interstorm period June 25, 1987 (day 176) to July 1 1 , 1987 (day 192). 


head value </f mjn to an arbitrarily small value (indicated as 
“-999.9" in Table 2) to ensure that evaporation remained 
atmosphere controlled throughout the simulation period, 
with the atmospheric evaporation rate equal to the actual 
(observed) rate. If potential evaporation data were available 
and an estimate of t/r min could be made, an important check 
on the performance of the model would be to run a simula- 
tion using the potential rates as input boundary conditions 
and allowing the model to switch from atmosphere- to 
soil-controlled evaporation, so that the computed actual 
evaporation rates could be compared with the observed 
actual rates. 

3.3. Simulation of 17-Day Rainfall-Interstorm Sequence 

The precipitation, evaporation, and streamflow observa- 
tion data used for the 17-day rainfall-interstorm simulation 
(June 25 to July 11, 1987) is shown in Figures 4 and 5. We 
used daily-averaged evaporation rates (Figure 5b) and had to 
interpolate the average rates for days 179 and 180 due to 
insufficient observation data for these two days. There were 
two major rain events during the 17-day period (Figure 4), on 
June 28 (17.14 mm) and July 5 (10.23 mm). In addition, 
between 2.3 and 3.5 mm of rain was recorded for three 
smaller precipitation events which occurred on June 29, June 
30, and July 7. Several storms also occurred during the week 
ending June 25, delivering 11.5, 21.2, 7.2, and 21.5 mm of 
rain on June 18, 20, 22, and 24 respectively. These large 
rainstorms immediately prior to June 25 enabled us to 
assume fairly wet initial soil conditions for the simulation. 



Time (day of year) 



Fig. 5. (a) Observed streamflow (solid line) and simulated discharge (dotted line), and ( b ) comparison of 

streamflow (solid line), evaporation (dotted line with crosses), and precipitation (dashed line) components of the 
catchment water balance for the 17-day rainfall-interstorm period ending July II, 1987 (day 192). 
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There were no observation data to verify the initial pressure 
head distributions used for the simulations, so the initial 
conditions can be considered somewhat arbitrary. We 
found, however, that varying the initial conditions did not 
have as great an effect on the simulated catchment discharge 
rates as variations in other model parameters. 

We ran two series of simulations using the 17-day rainfall- 
evaporation data. In the first series we used detailed, 15-min 
rainfall rates (Figure 4) and attempted to match the observed 
streamflow record. There are two distinct discharge events 
in the streamflow hydrograph (Figure 5 a ), in response to 
large storms on June 28 and July 5. Since the two discharge 
events appear to be responses to independent storms, wc 
calibrated the model to match the June 28 event and then 
used these calibration results to attempt a match of the July 
5 event. We remark that for proper model calibration and 
verification we should use additional observation sequences, 
encompassing a wide range of scenarios, though given the 
field data limitations described earlier we were unable to 
conduct more extensive tests. Furthermore, calibration 
based on streamflow observations alone is of limited use, as 
can be ascertained when we examine the relative contribu- 
tions of the various components of the catchment water 
balance. We see in Figure 5b that the rainfall and evapora- 
tion components are at least an order of magnitude larger 
than the streamflow component for the 17-day period of 
interest. It is apparent that almost all of the rainfall occurring 
during this period infiltrates the soil and only a small part 
reaches the stream as overland flow or rapid subsurface 
flow. 

The simulations of the 17-day period enabled us to exam- 
ine the sensitivity of model response to variations in param- 
eters and inputs. In particular, some preliminary tests on 
time aggregation effects were conducted in a second series of 
simulations. In these tests we used daily-averaged rainfall 
rates (Figure 5b), and we compared the hydrologic and 
numerical behavior with results from the simulations using 
detailed 15-min rain rates. We obtained significant discrep- 
ancies in surface saturation response during the storm 
events, but overall the simulated hydrologic responses using 
the averaged and detailed rain rates were quite similar. The 
advantage of using averaged rates for atmospheric inputs is 
that larger time steps can be used, making the simulation 
cheaper to run. Averaged rates also produce smoother 
boundary conditions which can reduce numerical difficulties. 
On the other hand, the high degree of variability in atmo- 
spheric fluxes may make it necessary to use detailed inputs, 
depending on the hydrologic responses of interest. For 
instance, large discharge events triggered by infiltration 
excess overland flow will be very sensitive to small-scale 
fluctuations in precipitation rates. 

The results from the first series of simulations are shown 
in Figure 5a. We were able to approximately match the 
observed peak, volume, and duration of streamflow re- 
sponse for the June 28-30 rainfall events, but the model 
produced no discharge for the July 5 storm. Note that the 
model result shown is total discharge, which is surface runoff 
routed to the stream. The model does not consider channel 
flow, and thus we could not compute streamflow rates at the 
catchment outlet. The lack of a channel flow component, 
along with the linear and discrete method used in the model 
to route and delay overland flow (surface runoff produced on 
30 x 30 m grid blocks is accumulated without downslope 


reinfiltration into 2-hour hydrograph intervals), may account 
for the fluctuations in the computed discharge shown in 
Figure 5a. The results for the first two days of the 17-day 
simulation are not shown in Figure 5a. There is streamflow 
activity on these first two days in response to the rain events 
which occurred during the week ending June 25, but we do 
not explicitly account for these rain events in the 17-day 
simulation. 

The model was first run using the parameter values 
obtained from observation data (Figure 3 and Tables 2 and 
3), but this produced an excessive amount of surface runoff. 
To calibrate the model we systematically altered the values 
of some of the parameters: soil depth L , vertical discretiza- 
tion Az, overland flow velocity, saturated conductivity K s at 
the surface, fitting parameter / for vertical exponential K s , 
time steps, and initial conditions. We found that only over- 
land flow velocity and surface K s had significant effects on 
the simulated discharge: overland flow velocity on the 
duration of the discharge response and surface K s on the 
peak rate and total volume of discharge. Vertical discretiza- 
tion will have some effect on simulated discharge in the case 
where infiltration excess runoff production is dominant. This 
type of runoff is controlled by surface K s , and in the model 
we evaluate relationship (5) at the midpoint of each element 
so that a thin surface layer will have a higher saturated 
conductivity than a thicker surface layer, thereby potentially 
reducing infiltration excess runoff. We found, however, that 
even using a 19-layer vertical discretization with the layer 
nearest to the surface 0.02 m thick (compared to 0.2 m for 
the four-layer discretization), the model still produced too 
much runoff. To reduce the amount of runoff it was neces- 
sary to increase surface K s from its fitted value of 0.0218 
m/h. We obtained best results using a value of 0.218 m/h (for 
the four-layer discretization). This higher value may be 
justified considering the remarks made earlier regarding the 
discrepancies between near-surface remotely sensed soil 
moisture observations and deeper measurements obtained 
from soil core analyses, since these discrepancies suggest 
that the soil near the catchment surface may be very porous. 
With surface K s = 0.218 m/h, an overland flow velocity of 
0.4 m/h, and all other parameter values as determined from 
observation data, we obtained the match shown in Figure 
5a. To match the streamflow response for July 5 it was 
necessary to run a separate simulation using surface K s = 
0.06 m/h. The need for different surface K s values to match 
separate discharge events suggests that it may be important 
to collect and incorporate information on the spatial distri- 
bution of hydraulic conductivity for the ID catchment. Note 
that since we altered surface K s from the value obtained by 
parameter fitting (Figure 3a) we should have refit parameter 
/ to the observation data using the modified value of surface 
K s . This was not done for the simulations reported here. 

The surface runoff which generated discharge during the 
17-day simulation was initially of the infiltration excess type, 
in response to peak rainfall rates on June 28, and later of the 
saturation excess variety, in response to declining rates but 
continued rainfall on June 29 and 30. This can be seen in 
Figure 6, where infiltration excess saturation is shown by 
heavily shaded areas, saturation excess by solid areas, and 
unsaturated portions of the catchment surface by lightly 
shaded areas. Note the high degree of surface saturation at 
the start of the simulation (day 176.0) owing to the wet initial 
conditions which were used. The initial conditions were 
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Fig. 6. Surface saturation response for the 17-day precipitation-evaporation simulation, using !5-min rainfall rates: 
(Top left) day 176.0 (June 25, 1987, initial conditions); (top right) day 179.103; (bottom left) day 181.280; and (bottom 
right) day 193.0. Light shading indicates unsaturated; heavy shading, infiltration excess; and solid regions, saturation 
excess. 


generated using a topographic index method as described 
previously, producing the wettest soils closest to the stream 
channel and catchment outlet. Since rainfall rates and satu- 
rated hydraulic conductivities were horizontally homoge- 
neous, infiltration excess runoff occurred first where near- 
surface pressure head gradients were lowest, corresponding 
roughly to the wetter near-stream regions. The patchiness in 
the infiltration excess pattern seen in Figure 6 is probably 
due to oscillations in the surface boundary condition switch- 
ing mechanism from one iteration to the next, oscillations 
which were caused by approximation errors in the back 
solved surface flux values. Saturation excess runoff occurred 
where water table levels were closest to the surface, corre- 
sponding once again to near-stream regions. These partial 
contributing areas then expanded upslope from the catch- 
ment outlet and outward from near-stream areas in response 
to continued rainfall. 

It is not surprising that we found high sensitivity to surface 


K s in the simulated discharge given that infiltration excess is 
a dominant mechanism for surface runoff generation. In fact, 
when we used daily-averaged rainfall rates we failed to 
generate any infiltration excess runoff, and consequently the 
simulated discharge hydrograph had a smaller peak and less 
volume and was of shorter duration than the observed 
streamflow hydrograph. The simulation with daily-averaged 
rates also produced less saturation excess overland flow. 
The different responses during periods of heavy rainfall can 
be seen in the vertical profiles of presssure head and mois- 
ture content shown in Figure 7. These vertical profiles are 
taken at the three selected surface nodes of Figure 2. In 
these plots we can see that the daily-averaged rainfall rates 
yielded drier soils during the storm on day 179. (The time 
value for the profiles labeled 179.1 in Figure 7 is actually 
179.105 for the 15-min rainfall rate case and 179.162 for the 
daily-averaged case; for the profiles labeled 188.5 the actual 
time is 188.461 for the 15-min case and 188.4% for the 
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daily-averaged case.) Aside from discrepancies in surface 
saturation results during heavy rainfall periods, the hydro- 
logic responses were generally quite similar between the 
simulations using averaged and detailed precipitation rates. 
This can be seen in Figure 7 for the profiles at day 188.5 and 
at the end of the simulation (day 193.0). 

The main difference in numerical behavior between the 
averaged and detailed simulations was that the simulation 
using 15-min rates required significantly more CPU time, 
since smaller time steps had to be used to resolve the 
detailed rainfall rates. The numerical results are shown in the 
convergence and mass balance plots of Figure 8 and are 
summarized in the second and third columns ffom the left in 
Table 4. The simulation using detailed rainfall rates had 
some difficulty during the peak rainfall/saturation periods on 
June 28 and July 5, as evidenced by the slower convergence 
and higher mass balance errors seen in Figure 8. Outside of 
these two periods the number of nonlinear (Picard) iterations 
remained roughly constant at two or three and the number of 
linear (conjugate gradient) iterations decreased steadily as 
the simulation progressed. The 17-day rainfall-interstorm 
simulation of the ID catchment required 19 min of CPU 
using 15-min rainfall rates and 6.5 min using daily-averaged 
rates. 

4. Sensitivity to Model Resolution 

4. 1 . Statement of the Problem 

In this section we examine the feasibility of using a 
detailed physically based model for catchment scale simula- 
tions. In particular, we will investigate computational issues 
related to model resolution by simulating the ID subcatch- 
ment over a range of grid sizes and convergence tolerance 
values. In a study of discretization clfccts for a two- 
dimensional hillslope model, Culver and Wood [1989J rec- 
ommended using A;r/Az ^ 20 for successful numerical 
simulations. In the results presented below we obtained 
satisfactory results using aspect ratios as high as 150 (A x - 
A y - 30 m, A z = 0.2 m). However, when we imposed a 
stricter convergence criterion on the nonlinear iterations, 
simulations with different aspect ratios produced different 
mass balance and convergence results. 

The ID subcatchment is small enough that it was not 
necessary to extrapolate the 30 x 30 m digital topographic 
data to coarser horizontal discretizations, but for larger 
catchments it will be necessary to examine whether aspect 
ratios even greater than 150 can be used. Large aspect ratios 
will also arise when a very fine vertical resolution is required 
to accurately reproduce runoff and moisture front responses 
from rainfall and evaporation events. Small vertical grid 
sizes may in turn lead to stability-related restrictions on time 
step, further increasing the cost of running large-scale catch- 
ment simulations. 

The effect of heterogeneities on numerical and physical 
grid constraints also needs to be investigated. For instance, 
while a fine grid may be needed to resolve hydrologic 
responses to highly variable inputs, a large problem size and 
high parameter variability may lead to ill-conditioned model 
equations which may adversely affect numerical conver- 
gence. Ababou et al. [1989] discuss related heterogeneity 
effects in the context of a three-dimensional finite difference 
mode! applied to steady state saturated flow problems. 


4.2. Simulation of 9-Day Interstorm Sequence 

In Figure 9 we show the evaporation and streamflow 
observation data for the 9-day interstorm period from May 
29 (day 149) to June 6, 1987 (day 157). There was no rainfall 
on the Kings Creek catchment during this period, except for 
a very small amount (I mm) on May 29. There was high 
streamflow just prior to day 149, and although we had no 
rainfall observation data prior to May 29 we assumed that 
there was significant rainfall immediately prior to this date, 
and we therefore started the 9-day simulation from saturated 
conditions. We used daily-averaged evaporation rates rather 
than detailed, diumally fluctuating rates. This allowed us to 
use large time steps (up to 12.0 hours), although we observe 
that the streamflow record in Figure 9 is oscillatory with a 
periodicity of approximately 1 day, which may be a response 
to the diurnal fluctuations in evaporation. Due to insufficient 
data for days 149, 153, and 154 we interpolated average 
evaporation rates for these days from the observed rates on 
surrounding days. 

We ran two series of simulations using the evaporation 
data. In the first series we simulated the entire 9-day period 
for a range of mesh discretizations from ID-30 to ID-03. This 
was done to study the effects of aspect ratio on catchment 
simulation responses, in particular to determine whether we 
can successfully run a catchment scale simulation using 
aspect ratios as large as 150 (corresponding to catchment 
ID-30). In the second series of tests we used only the first 12 
hours of the 9-day interstorm record and ran simulations of 
catchments ID-30 and ID-05 using a strict convergence 
criterion for the nonlinear iterations ( tol = 0.0005 m com- 
pared to tol = 0.05 m used for the 9-day simulations). 

The results from the first series of simulations are pre- 
sented in Figures 10 and 11 and summarized in the fourth, 
fifth, and six columns from the left in Tabic 4. In the figures 
we compare the ID-30 and ID-03 catchment simulations, 
and in the table we include also statistics from the ID- 10 
simulation. Similar results were obtained for catchments 
ID- 15, ID-06, and ID-05. 

In Figure 10 we show vertical profiles of pressure head and 
moisture content at the three selected surface nodes of 
Figure 2. There is reasonable agreement between the results 
obtained using a 30 x 30 m grid discretization (aspect ratio of 
150) and the results using a grid resolution of 3 x 3 m (aspect 
ratio of 15). We note also the similarities in the profiles at the 
three different locations owing to the spatial uniformity of 
atmospheric input rates, initial conditions, and parameter 
distributions. 

The differences in aspect ratio did not have significant 
effects on the numerical performance of the model. In the 
fourth, fifth, and sixth columns from the left in Table 4 we 
observe close agreement for the ID-30, ID- 10, and ID-03 
simulations. All three simulations were successfully com- 
pleted in 40 time steps, and none of the runs encountered 
difficulties severe enough to require back stepping. All 
simulations required between two and four nonlinear itera- 
tions to converge and an average of about 14 iterations to 
solve the linearized system of equations. In each case, 
convergence of the linear solver was slowest at the start of 
the simulation (20-25 iterations) and much more rapid by the 
end of the simulation (three to five iterations). Similar mass 
balance results were obtained for each of the simulations, 
with the largest discrepancies occurring during the first few 
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Time (day of year) Time (day of year) 

Fig. 8. Convergence and mass balance results for the 17-day precipitation-evaporation simulation, using 15-min 
(solid lines) and daily-averaged (dashed lines) rainfall rates. In Figure 8 d the peak value for the solid curve is -3.47%. 
(Day 176 is June 25, 1987.) 


time steps. These discrepancies at the start of the simulation, third time step (238% for catchment ID-03 and -40.5% for 

which can be seen in Figure 1 1, account in large part for the ID-30), values much larger than the errors computed at all 

differences in average mass balance errors reported in Table other time steps. We note that average mass balance errors 
4. High relative mass balance errors were obtained at the were significantly lower for the 17-day rainfall-interstorm 
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TABLE 4. Summary of ID Catchment Simulation Results 



17-Day Precipitation-Evaporation 
Simulation (Catchment ID-30) 

9-Day Evaporation Simulation 

12-Hour Evaporation 
Simulation 


15-min Averaged 
Rainfall Rates 

Daily-Averaged 
Rainfall Rates 

Catchment 

ID-30 

Catchment 

ID-10 

Catchment 

ID-03 

Catchment 
ID- 30 

Catchment 

ID-05 

Number of nodes 

1,570 

1,570 

1,570 

12,680 

135,355 

1,570 

49,295 

Number of time steps 

663 

223 

40 

40 

40 

81 

201 

Number of back steps 

1 

1 

0 

0 

0 

16 

50 

CPU, s 

1,150.3 

391.4 

66.84 

583.9 

6,452.6 

215.5 

19,597.0 

CPU/time step/node 

0.00111 

0.00112 

0.00106 

0.00115 

0.00119 

0.00169 

0.00198 

Nonlinear iterations/time 

2.26 

2.17 

2.40 

2.48 

2.53 

2.56 

2.64 

step 

Linear iterations/nonlinear 

11.08 

18.39 

13.92 

13.69 

14.78 

8.91 

7.47 

iteration 

Average absolute MBE|, m 3 

0.012 

0.031 

2.927 

3.044 

3.284 

0.00764 

0.00081 

Average relative MBE|, % 

0.038 

0.034 

2.021 

5.960 

7.255 

0.233 

0.070 


Convergence parameter tol = 0.0005 m. MBE denotes mass balance error. 


simulation than for the 9-day interstorm simulation (compare 
the second and third columns from the left in Table 4), owing 
to a smaller value of tol and smaller time steps used for the 
17-day simulation. 

The importance of not being limited by aspect ratio 
constraints in a catchment scale simulation is suggested by 
the storage requirements shown in Table 1 and by the CPU 
results shown in the fourth, fifth, and sixth columns, from 
the left in Table 4. With a grid resolution of 30 x 30 m we 
were able to simulate the 9-day evaporation sequence in 



Fig. 9. Evaporation (dotted line with crosses) and streamflow 
(solid line) observation data for the 9-day interstorm period May 29, 
1987 (day 149) to June 6, 1987 (day 157). 


slightly more than 1 min of computer time, requiring 0(1O 5 ) 
words of memory. The same 9-day simulation with grids of 
3 x 3 m required over 100 min of CPU and similarly 100 
times more memory. 

The results from the second series of simulations, using 
tol = 0.0005 m and the first 12 hours of the evaporation 
period, are summarized in the two rightmost columns in 
Table 4. These results can be roughly compared to the 
results in the fourth, fifth, and sixth columns from the left in 
Table 4 for the 9-day, tol = 0.05 m simulation. In reducing tol 
we have increased the accuracy of the simulations, decreas- 
ing by about 2 orders of magnitude the resulting mass 
balance errors. However, this gain is achieved at the ex- 
pense of a substantial increase in CPU. For catchment ID-30 
we simulated the 9-day evaporation period in 40 time steps 
and 67 s of CPU s with no occurrences of back stepping, 
while the 12-hour simulation required twice as many time 
steps, frequent back stepping, and over 200 s of CPU. For 
catchment ID-05 the 12-hour simulation required 5.5 hours 
of CPU. 

We observe from the two rightmost columns in Table 4 
that catchment ID-05 has more trouble converging than 
catchment ID-30 (201 time steps and 50 back stepping 
occurrences compared to 81 steps and 16 back steps for 
ID-30) and that it achieves smaller mass balance errors than 
ID-30. In this case the differences in aspect ratio (150 for 
catchment ID-30 and 25 for ID-05) had an effect on the 
convergence of the nonlinear iterations and on the accuracy 
of the simulation results (as reflected in the mass balance 
errors). 

5. Summary and Concluding Remarks 

We have described a three-dimensional physically based 
numerical model for the simulation of hydrologic processes 
at the subcatchment and catchment scales. The simulation 
model consists of a grid generator and a finite element code. 
The grid generator makes use of catchment information 
extracted from digital topographic data, and the numerical 
code is based on the nonlinear Richards equation for flow in 
variably saturated porous media. The model can be applied 
to catchments of arbitrary geometry and topography. Non- 
surface boundaries are considered impermeable, and grid 
spacing is uniform horizontally but can be variable along the 
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Fig. 1 1 . Mass balance results for the 9-day evaporation simulation, using 3 x 3m (solid lines) and 30 x 30 m (dashed 
lines) grid discretizations. (Hour 0 is May 29, 1987, 0000 hours.) 


vertical coordinate. The model automatically handles both 
soil-driven and atmosphere-driven inflows and outflows, and 
both saturation excess and infiltration excess runoff produc- 
tion. Atmospheric inputs can be spatially and/or temporally 
variable, and heterogeneities in hydraulic conductivity are 
also considered. Initial conditions can be generated from 
elevation and surface conductivity data, and extended van 
Genuchten equations are used to describe the nonlinear soil 
hydraulic characteristics. Simplifying assumptions are made 
to handle flow processes occurring on the catchment surface 
and in streams. Time stepping is adaptive and based on the 
convergence behavior of the nonlinear iterative scheme. 

The simulation model was applied to subcatchment ID of 
the Konza Prairie Research Natural Area in northeastern 
Kansas. Observation data collected during a 1987 field 
experiment were used to parameterize and calibrate the 
model. 

In simulations of a 17-day rainfall and interstorm sequence 
we attempted to match discharge produced by the model 
with observed streamflow. It was necessary to increase the 
fitted value of surface saturated conductivity to match the 
largest of two distinct streamflow events, but with this higher 
value of surface K s we were unable to match the smaller 
streamflow peak. Matching both discharge events simulta- 
neously may require information on the spatial variability of 
saturated conductivity for the ID catchment. The 17-day 
simulations were performed using both detailed 15-min rain- 
fall rates and daily-averaged rainfall rates. The detailed rates 
are necessary if we want to capture the surface saturation 
response of the catchment during periods of heavy rainfall, 


but in other respects the use of detailed and averaged rates 
produced similar hydrologic and numerical results. 

Simulations of a 9-day evaporation period were performed 
using a wide range of spatial grid discretizations in order to 
study model resolution effects. We obtained satisfactory 
results using a nonlinear convergence tolerance of 0.05 m 
and a grid aspect ratio as large as 150, indicating that 
horizontal grid dimensions may not be unreasonably con- 
strained by the typically much smaller vertical length scale 
of a catchment and by vertical discretization requirements. 
This is an encouraging result, although tests using heteroge- 
neous parameter distributions and even larger aspect ratios 
are recommended. Constraints on grid aspect ratios dictated 
by numerical stability, accuracy, and convergence require- 
ments also need to be investigated. 

Much of the work described in this paper is in preliminary 
stages. In future work we would like to extend the model to 
handle hysteresis and seepage faces, and improve on the 
simplifying assumptions used in runoff routing. For a more 
comprehensive model of catchment scale hydrologic pro- 
cesses, it will be necessary to couple the three-dimensional 
subsurface flow model to a physically based model of 
overland flow and channel flow, and to incorporate vegeta- 
tion, transpiration, macropore flow, and nonisothermal ef- 
fects. With more detailed and extensive observation data for 
the ID catchment (e.g., potential evaporation rates, surface 
soil moisture readings, and surface conductivity measure- 
ments) we will be able to run rigorous model assessment 
tests, including an evaluation of errors in the model formu- 
lation. 
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It would appear that catchment scale simulations can be 
feasibly performed using a detailed physically based model. 
The 9-day evaporation simulation of the 0.24-km 2 ID catch- 
ment required only 67 s of Cray Y-MP comuter time using a 
finite element grid of 1570 nodes, and 180 min of CPU using 
a 135,355-node mesh. From the fourth, fifth, and sixth 
columns from the left in Table 4 we estimate the CPU 
requirement for our simulation model to be 0.001 s (2.78 x 
10” 7 hours) per time step per node, using a relatively weak 
convergence tolerance ( to I - 0.05 m) and using the ITPACK 
routine SSORCG vectorized for the Cray Y-MP to solve the 
linearized system of equations. The U.S. Geological Survey 
digital elevation data for the entire King’s Creek catchment 
contain 12,913 pixels at 30 x 30 meter resolution. The 
surface area of this catchment is 1 1 .62 km 2 . A finite element 
discretization of the Kings Creek catchment using nine 
vertical layers would yield a 129,130-node mesh. A single 
100-time step simulation with this grid would require 3.6 
hours of CPU. The mesh for a horizontal discretization of 
15 x 15 m with four layers vertically would contain approx- 
imately 258,260 node, and would require just over 7 hours of 
CPU, again for a simulation of 100 time steps. Aside from the 
vectorized SSORCG solver, the numerical model was exe- 
cuted in scalar mode. Significant efficiency gains can be 
expected from vectorization and parallelization of other 
components of the model. One possibility would be to assign 
to each processor of a parallel computer the task of simulat- 
ing one subcatchment of a large catchment such as Kings 
Creek. 
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